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Abstract 


While solving contact problem, it is a common practice to use Hertz's expression 
for the contact area and a parabolic distribution of normal contact stress. How- 
ever, for problem involving large deformation of elasto-plastic contact bodies, these 
things are not applicable. A finite element code is developed for solving 3-D large 
deformation elasto-plastic contact problems. Frictional effects, however, are not 
included. A node-to-segment interface model is used for calculating the contact 
stiffness matrix Contact constraints are imposed by using the Lagrange multiplier 
method The material and geometric nonlinearities are handled by using updated 
Lagrangian method with a modified Newton-Raphson scheme. The code is vali- 
dated by solving some test problems. Finally a 3-D problem (of contact between 
a ball and plate) is solved for different combination of contacting materials to 
demonstrate the applicability of the method. 



Chapter 1 


INTRODUCTION 


Any mechanical load results from an interaction between two mechanical com- 
ponents, when they come into contact with each other. Contact interactions, 
thus, exist m virtually all structural and mechanical system. Some times these 
interactions are intentional like in metal forming, while sometimes they are 
accidental like in vehicle crash. 

Mechanical problems involving contacts are inherently non-linear. In large de- 
formation problems, additional non-linearities such as geometric and material 
nonlinearities may creep in. By nature, contacting surfaces involve friction 
phenomenon. Contact is usually a dynamics phenomenon but can be analysed 
as a static phenomenon if inertial forces are small. The actual contacting sur- 
face, the contact forces and the displacements of the contacting surface are all 
unknown prior to the solution of the problem. 


1.1 Literature Survey 


Study of contact problems began more than hundred years ago. Newton’s third 
law of motion and Coulomb’s friction law can be considered as early milestones 
in contact problems. With these one could solve rigid body contact problems 
in global sense. With the development of science of deformable bodies, local 
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phenomena like stress distributions on contacting surface started receiving 
attention. Hertz solved static contact problem in elasticity in the 1880s (Hertz 
1881,1882). He assumed contacting area to be small and contacting boundaries 
to be frictionless. Later on, Muskhehshvili (1975 ) and Gladwell (1973 ) solved 
some small deformation elastic contact problems using the analytical technique 
of integral equations. 

With development of computers, people resorted to numerical methods. Finite 
element method has been the most widely used method. In earliar studies, 
contacting bodies were assumed to consist of linearly elastic materials and 
hence were assumed to undergo only small deformation. The loading was 
assumed to be proportional. The node-to-node contact interface model (Fig. 
1 la) was used for discretising the contact bodies and the contact conditions 
( or constraints ) involving nodal displacements and forces were applied by 
modifying the combined stiffness or the flexibility matrix of the contacting 
bodies The investigators who used this model are Ohte (1973), Gartner (1977) 
(used proportional loading), Francavilla and Zienkiewicz (1979) (considered 
frictionless contact), Sachdeva and Ramakrishnan (1981a, 1981b) etc. 

Kalkar & Randen (1979), Hung & Sauxe (1980) and Mahmound et. al. (1991) 
used the technique of mathematical programming for solving the frictionless 
small deformaton elastic contact problem. They minimised the total strain 
energy to find the contact area and the contact pressure distributions. 

When the deformation is large, there is possibility that a pair of nodes, presently 
in contact, may undergo a large relative displacement. If node-to-node inter- 
face model is used, remeshing becomes necessary for carrying out the subse- 
quent analysis. To avoid frequent remeshing, one must use the node-to-segment 
interface model (figure 1.1b). When a node-to-segment interface model is used, 
the contact conditions acquires a complex form when expressed in terms of the 
nodal nodal varibles To impose these conditions, an additional set of finite el- 
ement equations, involving the contact stiffness matrix, needs to be developed 
using the principle of virtual work of the contact forces. The methods, which 
are used for implementing the contact constraints, are the Lagrange multi- 
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plier method, the penalty method and the transformation matrix. Kikuchi 
and Oden (1988) have provided the mathematical foundation for the Lagrange 
multiplier method and penalty method. 

(a) (b) 

Figure 1 1; (a) node-to-node interface model (b) node-to-segment interface 
model 

The work of Bathe & Chaudhary (1985) seems to be the first attempt to de- 
velop an expression for contact matrix using the above method. This was 
done for planer and axisymmetric problems assuming the contact segment to 
be straight. Lagrange multiplier method was used to impose the contact con- 
straints including Coulomb’s friction law. Non-linear contact kinematics was 
used for developing a consistent contact matrix first for 2-D problems by Wrig- 
gers k. Simo (1985) and then by Parish (1989). They used both the Lagrange 
multiplier method and penalty method. The transformation matrix method 
for 2-D problems was proposed by Chen and Yeh (1988). A formulation for 
non-linear frictional model proposed by Gallego & Anza (1989) using perturbed 
Lagrangian functional. All the above formulations are for static large deforma- 
tion elastic contact problem. Resently, some of these formulations have been 
extended to elasto-dynamic contact problem (Chaudhary and Bathe, 1990; 
Wriggers et al,1990) and dynamic elasto-plastic problems (Zhong, 1993) 



1.2 Objective of the Present Work 


The objective of the present work is to develop a simple yet reasonably ac- 
curate and efficient finite element formulation for a 3-D large-deformation 
elasto-plastic conatct problem without friction. Since, it is a large deformation 
problem, a node-to-segment interface model is used. A procedure described 
by Zhong(1993) is used to calculate the contact stiffness matrix because of 
its simplicity. Lagrange multiplier method is used to impose the contact con- 
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straints as penalty method is found to give reasonable solution for a limited 
range of penalty numbers. 

Elasto-plastic behaviour is modelled by a flow rule based on von-Mises yield 
condition and power-law type isotropic hardening Incremental analysis is 
carried out using the updated Lagrangian method (Bathe, 1996). Modifled 
Newton- Raphason method is used to solve the resulting finite element equa- 
tions Contact iterations are carried out within a Newton-Raphson iteration. 

A finite element code is developed which is first validated by solving some test 
problems Later a 3-D problem is solved to demonstrate the applicabilty of 
the formulation. 


1.3 Structure of the thesis 

The finite element formulation for Updated Lagrangian method and consti- 
tutive relationship between stress and strains for plasticity are developed in 
second chapter. The third chapter deals with calculation of contact forces and 
imposition of the contact constraints by Lagrange multiplier method. The 
fourth chapter deals with validation of computer code and solving ball over 
plate problem. The last chapter summarises the work carried out and suggests 
the possible further scope of the present work. 
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Chapter 2 


Finite Element Formulation For 
Elasto-Plastic Deformation 


2.1 Introduction 


When the deformation is large, displacements, strains and stresses, experienced 
by a body do not remain proportional to the applied load. The response of 
the body becomes nonlinear. The nonliearities involved in a large deformation 
elasto-plastic problem are of two kinds: geometric nonlinearity and material 
nonlinearity. For elasto-plastic materials, the Updated Lagrangian formulation 
is the most convenient as the stress-strain relationship for such materials is 
usually represented in an incremental manner. 

In this chapter, first the stress and strain measures appropriate for a large 
deformation problem are discussed. Later incremental equations for Updated 
Lagrangian method are derived. To deal with the material non-linearity, elasto- 
plastic model based on the isotropic hardening and the von-Mises yield criteria 
is developed. Finally, discretised finite element equations are derived. 
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2.2 Stress and strain measures for large defor- 
mation 


In a large deformation problem, configuration of the body changes continu- 
ously. We can’t calculate Cauchy stress at time t -I- At simply by adding the 
Cauchy stress at time t and the stress increment over the time At. Calcula- 
tion of Cauchy stress at time t + At must take into account the rigid body 
rotation of the material elements as the components of Cauchy stress change 
even with rigid body rotation. To deal with the above difficulty, we need to 
use appropriate stress and strain tensors, which do not change with rigid body 
rotation. Such tensors are called as objective tensors. 


One of the commonly used objective stress tensors in the Updated Lagrangian 
formulation is the second Piola-Kirchoff stress tensor (denoted by S^j). It is 
related to the Cauchy stres tensor (gij) by the relation 


t-f-At Q 




t+At. 


(t-f At ^i,m) ^ ^mn (t+At ^j,n) 


( 2 . 1 ) 


Here, t+At*^z,m represents the derivative of the position vector *5 at time t 
with respect to the one (‘"'■^*5) at time t + At. Further, ‘p and denote 
densities at time t and t + At respectively. The superscripts on S^, and cTij 
represents the deformed configuration at {t H- At) while the subscript on 
represents reference the configuration (at t). 

The corresponding work conjugate strain tensor is the Green-Lagrange strain 
tensor defined by 

^ ( 2 . 2 ) 

where ^ is the derivative of the displacement vector at time t + At with 

respect to the position vector at time t. 


In Eulerian formulation, where the rate of deformation tensor (symmetric part 
of the velocity gradient tensor) is used as the deformation measure, Jaumann 
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stress rate tensor (cTy) is used as the objective tensor. However, even in the 
Updated Lagrangian formulation, one can use At as the incremental ob- 
jective tensor. This tensor is related to Cauchy stress by the following 
relation. 

At — (T-ijAt o‘ii^(^dj}zAt'^ ^2.3) 

where the incremental rotation tensor is given by 

riyAt = — (Aujj — Au^^i). (2-4) 

and Atu represents the incremental displacement at time t 


2.3 Updated Lagrangian formulation 


In Lagrangian formulation, we follow the particles of the body in their motion, 
from the original to the final configuration of the body. One can find stresses 
and strains, either with reference to the original configuration, called the total 
Lagrangian formulation or with reference to the current configuration, called 
the Updated Lagrangian formulation. For an elasto-plastic problem, the Up- 
dated Lagrangian formulation is more convenient to use. This formulation 
[Bathe 1996], used in the present work is described in this section. 

Vitual work expression at time t+At for N bodies in contact can be expressed 
as follows 


E / (2.5) 

n=l 


where. 


n=l 

y [ 5 

^Jt+Ats, 


( 2 . 6 ) 
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Here, is the Cauchy stress tensor at time t+At, 5 is the virtual 

infinitesimal strain tensor at time t + At. is the virtual displacement 

at time t + At and is the domain at time t + At. is the part 

of the boundary on which the external forces act while is the contact 

boundary. The first term denotes the virtual work done by the applied 

forces while the second term is the virtual work done by the contact forces. 

For further derivations, the summation sign is omitted. The major difficulty 
in solving the above equation is that the configuration at time t + At is un- 
known. In the Updated Lagrangian formulation, the virtual work expression is 
transformed to an integral over the known domain at time t, using the second 
Piola-Kirchoff stress tensor and the Green-Lagrange strain tensor. 

Thus, the virtual work expresion becomes 

/ d^V = (2.7) 

The second Piola-Kirchoff stress tensor can be decomposed as 

+ ( 2 . 8 ) 

where, using equation 2.1, we can write 

ls„ = V., (2.9) 

Thus, 

= ( 2 . 10 ) 

Further, the Green-Lagrange strain tensor can be decomposed into linear and 

nonlinear increments as follows. 

5*+^\e,J=5Ate^, + SAtr)^, ( 2 . 11 ) 

where (from equation 2.2), 

At^ij = — (AtUij -F AtUj^i) (2.12) 

Atf]ij = — {AtUk^t) (AtU/.j) . 


(2.13) 
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Here, AtWjj stands for the derivative of Atu ( incremental displacement vector 
at time t ) with respect to tX (positon vector at time t ). Thus, equation (2.7) 
can be written with incremental decomposition as 

f AtS,,S{Ate,,)d^V + f AtS,,5{Atl^^J)d^V + [ V„<5(Atr7,,) 

J^v J^V 

+ d*V = (2.14) 

The 2nd integral on the left hand side is a higher order term compared to other 
terms and hence can be neglected. Using appropriate constitutive relation 
AiSij can be approximated as 

AfSij = tC'ijki Af6ki (2.15) 

Thus, equation 2.14 can be written as 

Ateki5{Ateij) d*V + J^^^alJ5{At'rj^J) <fV + 

V„(5(Aiey ) d^V = (2.16) 

Above equation can be employed to calculate the incremental displacement at 
time t + At, which then can be used to evaluate the incremental strain and 
stress. However, when these incremental displacement and stress are added 
to the their corresponding values at time t, it is observed that the total dis- 
placement and stress at time t At do not satisfy the virtual work equation 
(equation 2 5). That is, the internal virtual work (left hand side) dosen’t match 
with the external virtual work(right hand side). This out of balance virtual 
work occures as a result of linearisation performed. In order to reduce this out 
of balance virtual work, we need to perform an iteration in which the above 
solution steps are repeated until we reach a convergence. The equation 2.16 
then can be put in the iterative form as follows. In the A:-th iteration, the 
equation becomes 

(A>+4<4) '5(A.+A.et) + 

/ '+^VJ-^5(At+At4~^) (2.17) 
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The displacements in an iteration are updated as follows 

-t- Auf , (2.18) 

For the first iteration, the stresses, the strains and the constitutive tensors of 
the previous increment are used, thus 

— tCtjki (2.19) 

Equation 2.17 corresponds to Newton-Raphson iterations. In present analyses, 
the stiffness matrix (i.e. the first two terms on the left side) is not updated 
during the iterations, which helps in reducing the computational time. Such 
iterations are called modified Newton-Raphson iterations. 


2.4 Elasto-plastic model 


As stresses developed in a material exceed the yield stress, the linear rela- 
tionship between stress and strain no longer remains valid. We now develop 
relationship between stress and strain based on the von-Mises yield criteria 
and isotropic hardening. 


For an isotropic hardening material, plastic potential is given by, (Owen and 
Hinton, 1980) 

F{cr^j,p) = (TeqM - ay{p) ( 2 . 20 ) 


Note that 


F = 0 


( 2 . 21 ) 


represents the yield criterion. The plastic potential F depends on the Cauchy 
stress tensor a^j through it’s second invariant a^g called as equivalent stress 
and defined by 


^eq — 


I 2 

M/ 


( 2 . 22 ) 


where cr'j is the deviatoric part of cry. Further, F depends on the variable 
yield stress of the material, ay , through a hardening varible p. For the case of 
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strain hardening, p is identified as the equivalent plastic strain and hence 
defined as : 

p = e% = j (2.23) 

and 

= (2.24) 

Here, def^ is the plastic part of the incremental strain tensor de^j and the 
integration in equation 2.23 is to be carried along the particle path. The 
dependence of <7^ on p ( or ef 5 )is normally approximated by a power-law type 
of relationship 

<7, = K(el,T (2-26) 

Here, K is called the hardening coefficient and n is called as the hardening 
exponent. 


The incremental plastic strain (def^) is obtained from the plastic potential 
using the following relation : 

d<,=dA^ (2.26) 

where dA is a scalar. This equation is called as the flow rule. Differentiation 
of equation 2.20 with respect to gives 


^-_L ' 

dUij 2(Jeg^*'^ 

Using this one can determine dA as. 


(2.27) 


dA = def, 


(2.28) 


Further, the hardening 
used to express dA as : 


relationship and the yield condition (eq. 2.25) can be 


dA= ^ 
H 


dcr. 


eq 


H 


(2.29) 


TT _ 

^ A^P ■ 


where 


(2.30) 
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the slope of the hardening curve, can be obtained from equation 2.25. Sub- 
stitution of equations 2.27 and 2.30 in equation 2.26 leads to the following 
constitutive equation 


del = 


3 dcrgg / 


(2.31) 


This constitutive relationship between the deviatoric stress tensor and the 
incremental plastic strain tensor is not really convenient for the Updated La- 
grangian formulation for which the incremental stress-strain relationship is 
needed. This can be obtained from equation 2.31 as follows: 


def, = 


'-IJ 


2 dcTf^i 


Note that, from equations 2.20 and 2.26, we get 


da^_ ^ , 

dcTki dcTki 2a 


(2.32) 


(2.33) 


Substitution of equation 2.32 in equation 2.33 leads to the following incremen- 
tal plastic stress-strain relationship : 


j p ^ 


eq 


(2.34) 


The incremental elastic stress strain relationship is given by 

del = ^[-^dakkSij + (1 + ^^)da^^] (2.35) 

where de® is the elastic part of dcjj , E is the Young’s modulus and v is the 
Poisson’s ratio. Adding the two relationships, we get 


de,, = del + de^ 


1/ (1 -1- 1^) j. , , 

“I dikOji + 


> n 
tj^'kl 


Oa' (T: 




daki 


(2.36) 


This is the incremental elasto-plastic stress-strain relationship needed in the 
Updated Lagrangian formulation. However, it is the following inverse relation- 
ship which is more useful ; 


dcTij = Cfjndeki 


(2.37) 
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where 


riEP 


V 


^ l-2u 2(3/i + H)al 


eqj 


and ii is the shear modulus related to E and u by the relation: 

E 




2(1 + p) 


(2.38) 


(2.39) 


Note that, in equation 2.37, dcr^ cannot be interpreted as the incremental 
Cauchy stress tensor. Instead, it has to be identified as an increment of one of 
the objective stress tensors defined in section 2.2. 


2.5 Finite Element Formulation 


In this section, first the incremental elasto-plastic stress-strain relationship is 
put in a vector form Then, the incremental equation for the Updated La- 
grangian formulation is put in matrix a form. Finally, these equations are 
expressed in terms of nodal unknowns using suitable approximation for incre- 
mental displacement field. 


2.5.1 Incremental stress-strain relationship represent- 
ing material nonlinearity 


Since the stress-strain realtionship is path dependent for an elastic-plastic ma- 


terial, it is convenient to represnt it in incremental form, 
stress-strain relation is expressed in vector form as 

The incremental 

{da} = [C^^Kde} 

(2.40) 

where 


■(dcr} = {daxxj da^zi da^yj day^^ da^x} 

(2.41) 

and 


{d€}^ — {d^xxj dtyyi d(,ZZ1 ‘^dtxyt ‘idCyz^ ‘^d^zx} 

(2.42) 
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are respectively the vector forms of the incremental (objective) stress and strain 
tensors The expression for the elastic-plastic matrix [C^^] is derived from the 
plastic potential F. For its derivation, it is convenient to express the yield 
condition (eqs. 2.20 and 2.21) as 


- (^y(p) = 0 

(2.43) 

where 

{c} = {o'jca:) (^yyi ^xyi (^yz^zxi } 

(2.44) 

On the yield surface, dF = 0 and hence 


ydFrp dF 

(2.46) 

or 

{ay {da} -AdX = 0 

where, the flow vector {a} is deflned by 

(2.46) 

dF 

d{ay 

(2.47) 

the parameter A is given by 


= “dX a?'**’- 

(2.48) 


and dX is the same scalar which appears in equation 2.26. Equation 2.27 gives 
the following expression for the flow vector: 

(4 = ( 2 -«) 


where {cr}' is the vector form of the deviatoric part of cr,j. 

Equations 2.23 and 2.28 give 

dp = dX = de^g. (2.50) 

Substitution of equations 2.43 and 2.50 in equation 2.48 leads to the following 
expression for A: 

dCT y dcx y 


(2.51) 
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The total strain increment, after splitting into elastic and plastic parts, can be 
written in terms of the incremental stress vector by using the vector forms of 
the elastic stress-strain relationship (2.35) and flow rule(2.26) . Thus, 

dF 


{dc} = {de^} -f {de^} = [C]-\da} + dX 




= [C]-^{da} + dA{a} (2.52) 

Here, [C] is the matrix form of the elasticity tensor Ctjki • Premultiplying both 
sides of equation 2.52 by {a}^[(7] and using equation 2.46 , we get 

{ar[C]{de} 

A + {aFiCJia} 

Substituting the above expression for dX in equation 2.52 we get 


{de} = [C] W 


which leads to 


Therefore, 


{da} = [[C] - 
[Cf^=([C] 


A-F{aF[C]{a} 
[C]{a}{a}^[C] 


A + {a}T[C]{a}J 

[C]{a}{aF[C] 


{de} 


(2.53) 

(2.54) 

(2.55) 

(2.56) 


A + {a}T[C]{a}J 

where the definitions of o and A are given by 2.49 and 2.51. For a material 
with power strain hardening law (equation 2.25), the expression (2.51) for A 
reduces to 

A = Kn(^^r-^ (2.57) 


IS 


given by 










' l-u 

V 

V 

0 

0 

0 



V 

1-p 

V 

0 

0 

0 

[(7] = 

E 

V 

1/ 

\-v 

0 

0 

0 

(1 + z/)(l *— 2jy) 

0 

0 

0 

1-2!/ 

2 

0 

0 


0 

0 

0 

0 

\—2v 

2 

0 



0 

0 

0 

0 

0 

l-2i/ 

2 


(2.58) 


As stated before, {da} has to be interpreted as an increment of one of the 
objective stress tensors defined in section 2.2. 
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2.5.2 Finite Element Equations 

For developing the finite element equations, equation 2.17 need to be put in a 
matrix form. For notational convenience, equation 2.16 is used instead of 2.17 
for this purpose. Once finite element equations corresponding to equation 2.16 
are derived , the equation corrsponding to 2.17 can be obtained by carrying out 
the appropriate notational change. Owing to the symmetries in 
Atr)^J and cr^, equation 2.16 can be rewritten in the following matrix form: 

J.y ‘[r]{A,,} d‘V 

+ j (^{A,e}^) '{a} d‘V = ‘+“R. (2.59) 

where. 


t{Ae}- — A^Cyy, “^AfCxy, 2iAt6zx} (2.60) 

{o'} — { CTxxi ^yyi ^xyi ^yzi ^zx} (2.61) 

t{AT]} = AfU^y, AiU^zi AfV^xi AfV^y, AfV^z! ' ' '} (2.62) 

■‘[E] 0 0 

*[T] = 0 *[E] 0 (2.63) 

L 0 0 13. 

and 

^xx ^xy ^xz 

* [E] = *^y^ (2-64) 

^^zx ^zy ^ zz 


The matrix is the elastic- plastic matrix evaluated at time t which is 

given by equation 2.56. 

The domain is discretised into a number of elements and the incremental dis- 
placement field over a typical element e is approximated as 


t{Au} = [Q] t{Auy (2.65) 

where, 

t{AuY = {Atu^, Atv^, Atw^, Atu^, Atv^, }. 
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Here, ,AtW^ stand for the (unknown) incremental displacements of 

node in a; , y and z directions respectively and 



{QiF 1 

[Q] = 

{< 53 ^ J 


( 2 . 66 ) 


where 


{QiV = [ -^1 0 0 #2 0 0 $3 . . . ] , 

{ Q 2 f = [ 0 $1 0 0 ^>2 0 0 ...], 

{QzV =[ 0 0 $1 0 0 $2 0 ... ] 

(2.67) 


and $ 1 , the functions of ( x,y,z ), are called the shape functions. 


With displacement field defined by equation 2.65, the strain field within the 
element can be calculated in terms of nodal displacements. 


A.«., = ^ 


( 2 . 68 ) 


where 


t{Q^Z- ■ (2-69) 

Substituting equation 2.68 in equations 2.12 and 2.13 and arranging them in 
vector form, we get the strain displacement relations: 


,{A6} = ,|Bt] ,{A«}‘. (2.70) 

,{Ar,}=,[B„],{Au}' (2.71) 


where. 


ABl] = 


■ t{QiK 

t{Q2}5 

t{Qz}^y + t{Q2}% 
. + 1{<33}^ 


(2.72) 
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and 


,[BnY = { ,{Qi}l ,{«2K .W2}5 ■■■) (2.73) 

Using the strain displacement relations derived above, the contribution to the 
integral (eq. 2.59) from a typical element e (with volume ^Ve) can be written 
as 

t[C^^] tlBi] dV) + 

5t{AuY^ d‘y) t{AuY + (2.74) 

5t{AuY^ tlBif {a}* dV) = (t+At{F}^) 

Here, the contribution to the term *+‘^*72. from the element e has been expressed 
in terms of the elemental external force vectors t+At{FY using a standard 
procedure (Note that t+AtP^ includes both the applied as well as contact force). 
Since the variation in the incremental displacement vector is arbitrary 

and expressing the term within 1st parenthesis as linear elemental stiffness 
[ Kl Y f^hat of 2nd one as nonlinear elemental stiffness [ ]® the above 

equation can be written as 

{tlKiY + t[KNLY) 4- t{fY — t+At{-P’F (2.75) 

or 

t[KYt{^uY + t{fY = t+At{FY 

where 

t{/r= / 

is the internal force vector at time t. 

Assembling the elemental matrices [KY aud the elemental vectors t{fY 
t+At{FY over all the elements, we get the following global equation : 

t[K] t{Au} + ,{/} = t+At{F}. (2.78) 

Here t[F] is the global stiffness matrix, t{Au} is the global incremental dis- 
placement vector (at time t) and *{/} and t+At{F} are the global internal and 


(2.76) 

(2.77) 
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external force vectors respectively. Here, contribution to comes from 

both the applied forces as well as contact forces. Calculation of contact forces 
is discussed in next chapter. 


Equation 2.76 is the finite element version of equation 2.16. If, instead of 
equation 2.16 , equation 2.17 is used, the corresponding finite element equation, 
for A:-th iteration, becomes 

(2.79) 

In k-th. iteration, the stiffness matrix and the internal force vector 

are calcualted using the stresses of (A; — l)-th iteration and carrying 
out the integration over the updated volume of {k — l)-th iteration . Strictly 
speaking t+At{F} is also needs to be updated as the surfaces over which the 
applied and contact forces act also change. 

2.5.3 Solution Procedure 

As stated before, modified version of the Newton-Raphason method is used 
in this work. In this version, the stiffness matrix is not updated during the 
iterations. Thus, the solution procedure can be discussed as follows. 

In the current increment, first the stiflhiess matix is calculated based on the 
stress tensor obtained at the end of the the previous increment i.e. Vjj. This 
stiffness matrix is not changed till all the iterations of the current increment 
are over. For the first iteration of the current increment, is calculated 

using = Vij. Then equation 2.79 is solved, using the values of external 

and contact forces at time t + At to find {Au}F Then, incremental strain is 
calculated from equation 2.70 and incremental stress is determined using the 
relation 

,+i,{Aa}‘ = f [Cp'’{de}' (2.80) 

JAt 

Where the integration is carried out using Euler Forward integration[Bathe,1990]. 
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This is considered as the incremental Jaumann stress, a,. At. Prom this incre- 
mental Cauchy stress At is calculated using equation 2.3 and 2.4. This in- 
crement is added to to get . Next, the geometry and t+At{-P'} are up- 

dated using incremental displcement. Using the updated geometry and 
the internal force vector t+M{fV is calculated. Using t+At{F} - t+At{fV as 
the new right hand side, is obtained as the solution of equation 2.79. 

Now the incremental displacement becomes 

t{A«} = t{Au}^ + t{Au}2 (2.81) 

In the second iteration, internal force vector t+At{/}^ is calculated from 
using the procedure of previous paragraph . Further, t+Ai{-P'} 
is updated based on the geometric changes corresponding to t{Ait}^. Using 
t+At{-P’} — {fV as the right hand side, now t{Au}^ is determined by solving 
equation 2.79. The procedure is continued till the internal and external force 
vectors match within a tolerance, at which point iterations are stopped. This 
indicates that equilibrium has been achieved in the current increment. Then 
the calculation for next increment begins. 



Chapter 3 


Contact Formulation 


The quintessential part of a contact problem is to find the contact forces. 
The difficulty is, along the contact boundary, neither displacement nor force 
boundary conditions are specified. To find the contact forces, we have to 
impose certain contact constraints. 

In this chapter, first the contact constraints are discussed. Then, the ex- 
pressions for contact forces are derived using appropriate discretization of the 
contact boundary. Later, Lagrange multiplier method is discussed for im- 
posing the contact constraints. Finally, the algorithm for large deformation 
elasto-plastic contact problem is described. 


3.1 Contact constraints 

A typical two body contact system is shown in Fig 3.1. The bodies occupy 
the domains and ^ 'pj^g boundaries of 

and are denoted by *+^*5^ and and their interior volumes by by 

i+Atyi t+Ati [/2 respectively. 


central library 

1. 1. T., KANPUB 


4iitiKPh A 
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At any time At, the boundary of contact body n (n =1,2) can be partitioned 
as 


t+At (^n _ t+At y t+At qu y t+At gn 


(3.1) 


where, 


= boundary where displacements are specified 
= boundary where forces are specified 
= boundary where contact occurs 


Assuming the boundary *+■^*5'" to be smooth, outward unit normal vector 
at point on the boundary is denoted by Other two 

orthonormal boundary vectors are and . (See Fig. 3.2) 


Suppose that two boundary points and are in contact with each 

other at time t + At and the contact traction at is denoted as 

Then by Newton’s third law of motion, we have 


t+At^2 


(3.2) 


where can be expressed as 

3 

t+Atrj^n _ t-\rAt^n t+At jyn 
i=l 


(3.3) 
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Here, denote the conaponents of in the direction of 

If the two contact boundaries are not to be welded together, then the normal 
component of the contact traction cannot be tensile. Thus, we have 

t+At^n < Q 

This represents the traction constraint on the contact boundary. Besides the 
above condition, there are other types of constraints. Physical cosiderations 
require that particles of one body can’t penetrate the other body. Thus, for 
points in contact, penetration must be zero. This condition can be mathemat- 
ically stated as 

= (*+^‘52 _ . t+At^2 ^ Q ^3 5^ 

where is the penetration or gap between the points in contact. This 

condition is called as the kinematic constraint. 
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3.2 Contact Force Expressions 


As the contacting bodies are discretised for the purpose of developing finite 
element equations, the contact boundary gets automatically divided into sur- 
face elements, called as contact segments. While discretising the bodies, a 
node-to-segment interface model is used as the present work deals with a large 
deformation contact problem. Normally in a contact formulation, contact is 
considered only at discrete nodes. Here, it is assumed that a node of the hitting 
body (body 1) comes into contact with a segment of the target body (body 
2). The first step in the development of the contact force expression is the 
determination of the unit normal and tangent vectors at a point on the target 
segment. This is described in the next paragraph. 

Geometry of a target segment can be described using 2-D shape functions. 
Thus, 

>)) = 'E K , (3.6) 

where (^, rj) are the shape functions and denote the nodal position 

vectors. 

The side of the segment on which contact can occur is refered to as positive side 
of the segment and the other side as negative side of the segment. The local 
numbering of nodes on a contact segment should always be counter clockwise 
when we face the positive side of the segment. This is to get the correct sense 
for the outward normal vector at a given point. A normal vector at a point 
*+^*5 (^, 77 ) on the positive side of segment is defined as 



(3.7) 

where, and are two tangent vectors at 

f+Af,j2(^ ,,) giyen by, 

<1 

II 

<I 

4- 

(3.8) 

'+“JVJ=*+A‘x2K,,)_, 

(3.9) 

The unit normal and tangent vectors are defined as 

t+Atjy2 _ t+At^2 

(3.10) 
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«+At^2 ^ t+Atjy2/|i+Atjy2| 

*+^*iV| = (3.12) 

Since contact is considered only at discrete hitting nodes, we deal only with 
point forces at discrete hitting nodes rather than with distributed contact 
tractions. At time t+At, the ]|^g array of the cartisian components 

of the point force vector at a point of the target segment at which the 

hitting node makes the contact. By Newton’s third law, the force vector on 
the hitting node body will be —p+^*F^}. Let and be the 

virtual displacement vectors (at time t + At) at these points ( on the target 
and hitting bodies respectively ) when they come in contact. Then the virtual 
work at hitting node can be expressed as 

5wc = ^t+Atp2y 

Note that the elements of namely are the 

(virtual) nodal displacements of the hitting node. The elements of 
can be expressed in terms of the (virtual) nodal displacements of the target 
segment z = 1, - • ■ N) using 2-D shape functions 

{(pi t = l,N). Thus 

_ ^t+At^ly ^ [Q,]{S^+^*Uc} ( 3 . 14 ) 

where 

—1 0 0 * 1*1 0 0 $2 0 0 $iv 0 0 

Qc = 0 -1 0 0 #1 0 0 $2 0 • • • 0 0 

[ 0 0 -1 0 0 #1 0 0 <^2 0 0 

(3.15) 

and 


St+At^l 


(3.16) 
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The vector can be expressed as 



- t+At^2^ ■ 



> 

to 

II 

t+Ai^2^ t+At^2^ t+Atj^2^ 


‘+'"71 - 


t+AtN2^ ‘+AtJV|3 *+AtiV|3 


t+At f 2 

1 /3 J 




Where represents the j-th cartesian component of the unit vector 

and are the components, with respect to of the con- 

tact force vector *+^*#2 Combining equations 3.13, 3.14 and 3.18 we get the 
following expressions for the virtual work at the hitting node: 

Swc = (3.19) 

The vector at the hitting node is given by 

{‘+^Vc} = [QcY [*+A*iV2j|i+Aty2| (3 20) 

The vector can be split as a contribution from normal contact forces 

and tangential contact forces. Thus 


Where 


{‘+AVj = (*+AVj + {*+AV^} 

^ [Qjr ^t+Atj^2y+Atf2 


(3.21) 

(3.22) 


and 


{‘+AV_^} = [Q,f ({*+^*iV|}*+At_^2 ^ |t+Ai^2^t+A*_^2) (3 33) 

Here the arrrays z = 2,3 contains the cartesian components of the 

vector *+AZ]y 2 _ 


Before discussing the method of determining the contact forces, it is necessary 
to express the kinematic constraint in terms of the nodal displacements. The 
matrix form of equation 3.5 is 

*+Atp ^ {t+At pj2yT y+At^2 _ t+At^l y (3 24) 
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Note that 


^t+At^2 _ ^ (3 25 ) 

where the second term on the right hand side is an array of the Cartesian 
componenets of the difference of the incremental displacement vectors of the 
target and hitting bodies at the contact node. Analogous to equation 3.14 this 
can be written as 

{Atu^ - Atu^} = [Q^{AtUc} (3.26) 

where the vector 

{AtUc}'^ = {Atu\ Atv\ Atw^, Atul, Atvl, Atw\, • • • , twh} (3.27) 

contains the nodal values of the incremental displacement of the hitting node 
and the nodes of the target segment. Recognizing the first term of the right 
hand side of equation 3.25 as (penetration at time t), equation 3.25 now can 
be written as 

= V + [Qc]{tAUe} (3-28) 


3.3 Lagrange Multiplier Method 


In Lagrange’s method, the contact forces are considered as primary unknowns 
and the kinematic contact constraint is enforced exactly. In the present work, 
the contact problem is considered as frictionless. Thus contact force has only 
one component, namely, normal force component. Now equation 3.20 becomes 

{*+AVe} = [QcY (3-29) 


The kinematic constraint corresponding to the normal contact force is that the 
penetration should be zero. Thus, from equation 3.28, we have 


{‘+“jv?y[Qj{A,«„}+‘p=o 


( 3 . 30 ) 
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Combining equations 3.29 and 3 30, we obtain 


0 ■ 


[ {AtuJ 1 


r {‘+“re} 1 

_ 0 

< 

1 1 

i J 

^ =: ^ 

1 -v J 


(3.31) 


where 




(3.32) 


For L contacting nodes, there will be L such equations. Assembling these 
equations, we get the following global equation. 


0 




{AtC/J 

^t+Atp2y 


. -m 


(3.33) 


This equation is combined with equation 3.79 of chapter 2 to find the nodal 
displacement vector (which contains the vector {AtUc} of the contact displace- 
ments) and the contact force vector }. Note that the vector 

is a part of the vector t+At{F} on the right hand side of equation. 


3.4 Algorithm for elasto-plastic contact prob- 
lem 

For solving a typical elasto-plastic contact problem, the following steps are 
used. 

1. Contact searching: Contact node searching is carried out to find the 
potential nodes that (for both hitting and target bodies) may come in 
contact with the application of incremental load. The contact searching 
is carried out using a master and slave algorithm (Zhong 1993). 

2. Renumbering of nodes: Nodes for the hitting and the target bodies are 
renumbered such that contact nodes are numbered first. This is carried 
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out to facilitate static condensation of stiffness matrix which helps in 
reducing computational time. 

3. Stiffness matrix: Based on the current geomentry and state of stress 
stiffness matrix is calculated. 

4. Newton Raphason iteration: Newton-Raphason iteration counter is ini- 
tialised as i = 1. 

5 Initial penetration: Initial penetration *+'^*p* ^t+AtpO _ tp'j jg found using 
current geometry. 

6. Contact iterations: Within Newton-Raphason iteration contact itera- 

tions are carried out to find the actual nodes in contact. The constrant 
t+Atf 2 ^ Q is similar to equation 3.4) is used to determine whether 

the assumed node is in contact or not. The, displacements and contact 
forces are calculated after the convergence of contact iterations. 

7. Updating: Displacements, stresses and the contact forces are updated. 

8. Convergence: Newton-Raphason iterations are checked for convergence 
as described in chapter-2. If convergence is achieved then equilibrium 
state is reached else we have to go to step 5 and iterate till we get 


convergence. 



Chapter 4 


Results and Discussion 


A 3-D finite element code is developed for a large deformation elasto-plastic 
frictionless contact problem based on the formulation given in chapter-2 and 
chapter-3. 8-noded brick elements are used for discribing the contacting bodies. 
8-point Gauss Legendre numerical scheme is used for evaluating the stiffness 
matrix. Convergence tolerance used is 1 % on the global force vector. First 
some test problems are solved for validating the code. Later a problem of ball 
over a plate is solved, to demonstrate the applicabilty of the methos. 


4.1 Validation of computer code: 

Validation of code is carried out in two steps 

1. Validation for a large deformation elasto-plastic single body problem. 

2. Validation for a large deformation elastic frictionless contact problem. 
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4.1.1 Validation for elasto-plastic part of computer code 


The elasto-plastic part of the computer code has been validated with a theo- 
ritical example given in Chakrabarty (1987). The problem consist of a simply 
supported beam with uniformaly distributed load(figure 4.1). 



Figure 4 1; Simply supported beam under uniformlly distributed load 

Plane strain condition is simulated by applying appropriate boundary condi- 
tions. Only a half of the beam is analysed. 

Stress-strain relation used is 

{ojoy) = {Ee^lcjyY for a > Gy (4 1) 


Where the material constants are: 

Yield stress {cy) = 21.5 N/mm2 
Young’s modulas {E) = 2500 N/mm2 
Hardening coefficient (n) = 0.2 

Dimensions of the beam solved are 

Length (L) = 100.0 mm Breadth (B) = 10.0 mm Height (H) = 10.0 mm 

The ratio of maximum load to load at yield is plotted against the ratio of 
maximum deflection to deflection at yield (figure 4.2). The results are in good 
agreement with theoritical solution provided by Chakrabbrty (figure 4.3). 
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4.1.2 Validation of contact part of code 


The contact part of the code has been tested with a numerical problem solved 
in Zhong (1993). The problem consists of two elastic cantilever beams in 
frictionless contact (figure 4.4). 



i 

H 

h 

H 

T 

MN) 


Where dimensions are 
length (L) = 1.0 m 
height (H) = 0.1 m 
length (B) = 0.1 m 
gap (h) = 0 01 m 
Material properties used are: 

Young’s modulus {E) = 200 GPa 
Poissonis ratio = 0.0 

Zhong (1993) has solved it as a 2-D plane stress problem. In present analysis 
plane stress condition is simulated by applying appropriate boundary condi- 
tions. The finite element mesh used is shown in figure 4.5. Two elements have 
been chosen in the direction perpendicular the paper. 

The Lagrange multiplier method is used to impose contact constraints. The 
total load is applied in thirty steps. The two deformed configurations are 
shown in the figure 4.6 and 4.7. 
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Figure 4.6: Deformed configuration of the beams in contact (Load level-15) 



Figure 4.7: Deformed configuration of the beams in contact (Load level-30) 
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The displacement of the node at which load is applied is plotted against the 
load level in figure 4.8. The results are in good agreement with the results 
provided by Zhong (figure 4.9). 



Load level 


Figure 4 8: The maximum deflection of the upper beam as a function of the 
load level 



Figure 4.9: The maximum deflection of the upper beam as a function of the 
load level 
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4.2 Elasto-plastic contact problem 


In this section, an elasto-plastic contact problem between a ball and a plate is 
solved. Geometry of the contact system is shown in figure 4.10. 



//TTTTTV 






Figure 4.10; Contact system of a ball over a plate (R = 75.00mm, h-20.00mm) 


By exploiting the symmety in two perpendicular planes, only a quarter of the 
domain is analysed. The finite element mesh used, in a plane of symmetry is 
shown in figure 4.11. For the ball five elements are chosen in the Q direction 
while for the plate, eight elements are selected in the third direction. 

The ball and plate are made of the polycarbonate material. The material 
properties are: 

E (Young’s modulus) = 1114.58 N/mm^ 

V (Poisson’s ratio) = 0.36 
Uy (Yiled stress) = 44.58 N/mm^ 
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Figure 4.11: Finite element mesh for the ball over plate system 
No of nodes for body ball = 100 

af { Failure stress) = 82.71 N/mm^ 
k (Hardening coefficient) = 87.72 N/mm^ 
n (Hardening exponent) = 0.909 


Analysis is carried out till the plate fails. The deformed configuration at failure 
is shown in figure 4.12. 

Initialy, only one row of hitting nodes is in contact with the plate. As the load 
is increased, more and more nodes of the hitting body come in contact with 
the plate. The maximum deflection occures at the point of application of load 
while the element just below the applied load is maximally stressed. Figure 
4.13 shows the variation of load level versus the maximum deflection. 

To study how the contact length varies with different materials, following com- 
binations of materials are used for the ball and the plate. 


1. Acrylic ball and acrylic plate 
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Figure 4.12: Deformed mesh (polycarbonate-polycarbonate) 



Load level 


Figure 4.13: Maximum deflection of load point as function of load level 
(polycarbonates-polycarbonate) 


2. Acrylic ball and polycarbonate plate 

3. Polycarbonate ball and acrylic plate 


The material properties of the acrylic are: 
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E (Young’s modulus) = 2183.99 N /mm^ 
p (Poisson’s ratio) = 0.35 
ay (Yiled stress) =21.5 N/mm^ 
af { Failure stress) = 73.39 N/mm^ 
k (Hardening coefficient) = 1408.6 N/mm^ 
n (Hardening exponent) = 1.01 


Figures 4.14, 4 15 and 4.16 show the deformed configurations at failure for dif- 
ferent material combinations. From these figures, one can see that in the case of 
acrylic-acrylic combination only one row of hitting nodes come in contact while 
for the other two combinations, two rows of hitting nodes come in contact. 
Figures 4.17, 4.18 and 4.19 are the plots of load level against maximum the 
deflection. It is found that in case of acrylic-acrylic and polycabonate-acrylic 
combinations plate fails before the ball while in polycarbonate-polycabonate 
and acrylic-polycarbonate cases the ball fails before the plate. From the load 
deflection plots, one can see that curve is more stiffening for the case of 
polycarbonate-polycarbonate (4.13) and acrylic-polycarbonate (4.18) combi- 
nation. This is due to the fact that though the deflection is more in these 
cases, the plastic deformation zones are smaller than the other two cases. This 
may be attributed to the fact that the yield stress of polycabonate is large 
than of acrylic, hence many points do not yield. 

One of the important objectives of an analysis of a contact problem is to de- 
termine the contact stress distribution. A fine mesh is required in the contact 
region for accurate estimation of the contact stresses. In this work, the refine- 
ment of mesh in the contact region can not be carried out due to difficulties 
involved in generation of a large 3-D mesh. 
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Chapter 5 


Conclusions and Scope for 
Future Work 


5.1 Conclusions 


A finite element code has been developed for solving 3-d large deformation 
elasto-plastic frictionless contact problems. A node-to-segment inerface model 
is used for developing the contact stiffness matrix. The Langrage multiplier 
method is used for inforcing the contact constraints. The updated Langrangian 
formulation along with modified Newton-Raphson iterative scheme is used for 
carrying out the incremental analysis of the elasto-plastic deformation. 

The code has been validated by solving a few test problems. Finally, a 3-D 
problem (of contact between a ball and plate) is solved for different combina- 
tions of contacting materials. The result seems reasonable. This shows that 
the method can be used successfully for solving 3-d large deformation elasto- 
plastic (frictionless) contact problems. 


5.2 Scope for Future Work 


The present can be extended in the following directions : 
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• Friction effects can be included. This involves modification of the contact 
stiffness matrix. 

• To include the code for solving impact problems, dynamic (inertial) ef- 
fects should be included. A finite difference scheme needs to be incorpo- 
rated to tackle the dynamic terms. 

• An effective pre-processor should be added to the code so that the prob- 
lems with complicated geometry can be solved, a versatile post-processor 
is also desirable to represent the displacement and strss fields in the con- 
tact region. 
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